Beating the efficiency of both quantum and classical simulations with semiclassics 
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While rigorous quantum dynamical simulations of many-body systems are extremely difficult (or impossible) 
due to the exponential scaling with dimensionality, corresponding classical simulations completely ignore quan- 
tum effects. Semiclassical methods are generally more efficient but less accurate than quantum methods, and 
more accurate but less efficient than classical methods. We find a remarkable exception to this rule by showing 
that a semiclassical method can be both more accurate and faster than a classical simulation. Specifically, we 
prove that for the semiclassical dephasing representation the number of trajectories needed to simulate quantum 
fidelity is independent of dimensionality and also that this semiclassical method is even faster than the most ef- 
ficient corresponding classical algorithm. Analytical results are confirmed with simulations of quantum fidelity 
in up to 100 dimensions with 2 1700 -dimensional Hilbert space. 

PACS numbers: 05.45.Mt, 03.65.Sq, 05.45.Pq, 05.45.Jn 



Introduction. Conect description of many microscopic dy- 
namical phenomena, such as ultrafast time-resolved spectra or 
tunneling rate constants, requires an accurate quantum (QM) 
simulation. While classical (CL) molecular dynamics simula- 
tions are feasible for millions of atoms, solution of the time- 
dependent Schrodinger equation scales exponentially with the 
number D of degrees of freedom (DOF) and is feasible for 
only a few continuous DOF. An apparently promising solu- 
tion is provided by semiclassical (SC) methods, which use CL 
trajectories, but attach to them phase information, and thus 
can approximately describe interference and other QM ef- 
fects completely missed in CL simulations. Unfortunately, SC 
methods suffer from the "dynamical sign problem" due to the 
addition of rapidly oscillating terms, resulting in the require- 
ment of a huge number of CL trajectories for convergence. 
Consequently, most SC methods are much less efficient than 
CL simulations and in practice were used for at most tens of 
DOF. Even though several techniques have explored this is- 
sue Oj], the challenge remains open. Below we turn this chal- 
lenge around by showing that in simulations of QM fidelity 
(QF) Bi, a SC method called "dephasing representation" 
(DR) is not only more accurate but, remarkably, also faster 
than the most efficient conesponding CL algorithm 0. 

Quantum and classical fidelity. QF was introduced by 
Peres |5|] to measure the stability of QM dynamics (QD). He 
defined QF Fqyiif) as the squared overlap at time t of two QM 
states, identical at t = 0, but subsequently evolved with two 
different Hamiltonians, Hq and H e = Hq + eV: 



FQu(t) := |/qm(*)| i 
/qm(£) := {^U-'U'M), 



(1) 
(2) 



where /qm(£) is the fidelity amplitude and XJ\ := 
cxp(—iH e t/ti) the QM evolution operator. Rewriting Eq. (O 
as /qm(*) = (i> \U l | tp) with the echo operator U l := U^Uq, 
it can be interpreted as the Loschmidt echo, i.e., an overlap of 
an initial state with a state evolved for time t with Hq and 
subsequently for time — t with H t . (In general, we write time 
t as a superscript. Subscript e denotes that H e was used for 



dynamics. If an evolution operator, phase space coordinate, 
or density lacks a subscript e, Loschmidt echo dynamics is 
implied.) QF amplitude (O is ubiquitous in applications: it 
appears in NMR spin echo experiments [6], neutron scatter- 
ing 01, ultrafast electronic spectroscopy jM l^l, et c. QF (Q} is 
relevant in QM computation and decoherence itlQfl . and can be 
used to measure nonadiabaticity 111 10 or accuracy of molecular 
QD on an approximate potential energy surface III 211 . 

Definition (Q]i can be generalized to mixed states in different 
ways I2L 1 1 314151 . but we assume that the initial states are pure. 
In this case, one may write QF (HJ as Fq M (£) = Tr (/3'p ) 
where p* := U^pU^ 1 is the density operator at time t. In 
the phase-space formulation of QM mechanics, QF becomes 
F QM (t) = h~ D J dxpl w (x)p t o w (x) where x := (q,p) is a 



point in phase space and Aw(x) :— J d£,{q — £/2 



.4 



suggests its CL limit, called CL fidelity (CF) IUrjl[l7ll 



is the Wigner transform of A. This form of QF 



FcUt) ■= F M (t) 



dxp\{x)pl(x) 
dxp t (x)p°(x) 



(3) 



(4) 



where the first and second line express CF in the fidelity and 
Loschmidt echo pictures, respectively. If F or p lack the sub- 
script "CL", "QM", or "DR", "CL" is implied. 

Dephasing representation. There were several attempts at 
describing QF semiclassically. Most were analytical yl [Till 
and valid only under special circumstances because the nu- 
merical approaches were overwhelmed with the sign prob- 
lem. Extending a numerical SC method for localized Gaus- 
sian wavepackets (GWPs) 11911 . the DR was introduced as a 
more accurate and general approximation of QF lfT3Ul5ll . The 
DR of QF amplitude is an interference integral 



A*(t) := h~ D J <fcV(s°)exp[t0(a: o ,t)], (5) 
b(x°,t) := -AS(x°, t)/h = (e/h) [ &rV(x T e/2 ), (6) 



2 



where the phase is determined by the action AS due to the 
perturbation along a trajectory propagated with the average 
Hamiltonian H t/2 HHq]. Above, x\ := $'(a; ) where 
is the Hamiltonian flow of and the perturbation V can, in 
general, depend on both q and p. The DR of fidelity, com- 
puted as For := |/dr| 2 , was successfully used to describe 
stability of QD in integrable, mixed, and chaotic systems 11131 - 
H5I1 . nonadiabaticity pit] and accuracy of molecular QD on an 
approximate potential energy surface lfl2ll . and the local den- 
sity of states and the transition from the Fermi-Golden-Rule 
(FGR) to the Lyapunov regime of QF decay [21]. The same 
approximation was independently derived and used in elec- 
tronic spectroscopy (8|] . Recently, the range of validity of the 
DR was extended with a SC prefactor 02011 . The remarkable 
efficiency of the original DR observed empirically in appli- 
cations led us to analyze this property rigorously here and to 
compare it with the efficiencies of the QM and CL calcula- 
tions of QF. 

Algorithms. The most general and straightforward way to 
evaluate Eqs. (O-© and (0 is with trajectory-based methods. 
While the DR (0 is already in a suitable form, Eqs. ©-(lUl for 
CF must be rewritten using the Liouville theorem as 



F M (t) 



dx°p(x e t )p{x *) and 



echo V 



(t) = br D / dx°p(x- t )p(x°). 



(7) 



(8) 



Above, x l := where $* := $7* o $^ is the Loschmidt 

echo flow. Since it is the phase space points rather than 
the densities that evolve in expressions ©-(UK we can take 
p = Pyf°. For numerical computations, Eqs. © and (0)-<[8j 
are further rewritten in a form suitable for Monte Carlo eval- 
uation, i.e., as an average 

_ fdx°A(x°,t)W(x ) 



(A(x°,t)) 



W(x°) 



fdx°W(x°) 

where W is the sampling weight for initial conditions x°. Us 
ing W — pw(x°), the DR algorithm becomes 113141511 

/ DR (i) = (exp [i<j)(x°,t)]) pvi{x0) . 



(9) 



Sampling is straightforward for /?w > 0, but can be done also 
for general pure states B15I1 . Whi le p reviously used CL al- 
gorithms sampled from W = p &lH I22I1 . Ref. J3l consid- 
ered more general weights W — Wm(x°) := p(x°) M and 
W = W M (xo *) = p($(7*(a; )) M for the echo and fidelity 
dynamics, respectively. These weights yield four families of 
M-dependent algorithms |4] 



F M . M {t) 

Fecho-M(t) 
-Ffid-N-Af (*) 

N-m(*) 



t\l-M\ 



{p{x 7 >)p{xjy- M ) p(x - t) 
<p(0 2 - M > 

(p(ar-*)p(aj ) 1 - M > 



p{x°) M 



2-M\ 



(10) 

(11) 

(12) 
(13) 



where Im '■= h~ D J p(x°) M dx° is a normalization factor. 
Conveniently, the "normalized" (N) algorithms (fT2b-(fT3b do 
not require the normalization factor Im which is, for general 
states, known explicitly only for M 6 {0, 1,2} (I = nf , 
I\ = I2 = 1)- For further details, see Ref. J4J] where it was 
found that the echo-2 algorithm is optimal since it is already 
normalized (i.e., echo-2 = echo-N-2), applies to any pure state 
(in particular, p does not have to be positive definite), and- 
most importantly-is by far the most efficient CL algorithm. 

Efficiency. The reader does not have to be persuaded of the 
exponential scaling of QD with D. We just note that the direct 
diagonalization of the Hamiltonian leads to a QD algorithm 
with a cost O(t n 3 D ) — O{t n\ D ) where nn — n,f is the di- 
mension of the Hilbert space of D DOE Despite the indepen- 
dence of t, the scaling with D is overwhelming. More practi- 
cal is the split-operator algorithm requiring the fast Fourier 
transform (FFT) at each step. The complexity of FFT is 
0(riD logrio), hence the overall cost is 0(iDn\ D logrti). 
The effective n\ is reduced in increasingly popular methods 
with evolving bases, but the exponential scaling remains. 

Regarding the CF and DR algorithms, efficiency of 
trajectory-based methods depends on two ingredients: First, 
what is the cost of propagating N trajectories for time tl Sec- 
ond, what TV is needed to converge the result to within a de- 
sired discretization error crdjacr? As this analysis was done for 
the CL algorithms in Ref. [4], here we only outline the main 
ideas and apply them to analyze the efficiency of the DR. 

The cost of a typical method propagating N trajectories for 
time t is 0(c[tN) where Cf is the cost of a single force evalu- 
ation. However, among the above mentioned algorithms, this 
is only true for the fidelity algorithms with M = (i.e., fid-0 
and fid-N-0) and for the DR! Remarkably, in all other cases, 
the cost is 0(c(t 2 N). The cost is linear in time for a single 
time t, but becomes quadratic if one wants to know CF for all 
times up to t. For the echo algorithms, it is due to the neces- 
sary full backward propagation for each time between and 
t. For the fidelity algorithms, it is because the weight func- 
tion p(x^ t ) M changes with time and the sampling has to be 
redone for each time between and t |4J] . 

The number N of trajectories required for convergence can 
depend on D, t, dynamics, initial state, and method. Below 
we estimate N for the DR analytically using the technique 
proposed in Ref. |4]. The expected systematic component 
of (Tdiscr is zero for /dr and 0(iV _1 ) for For and is negli- 
gible to the expected statistical component a = 0{N^ 1 / 2 ) 
which therefore determines convergence. Expected statistical 

error of A(N) is computed as a\ (N) = \A(N)\ 2 -\A(N) 2 
where the overline denotes an average over infinitely many 
independent simulations with N trajectories. 

The discretized form of Eq. (0 is foR(t,N) = 

from which |/dr(£, iV)| 2 = 
2 

N- 1 



+ (l-N- 1 )F DR (t), 



F DR (t), and aj D 



A r_1 (l — Fur). The analogous calculation for Fur is some- 
what more involved but straightforward. Inverting the results 
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FIG. 1. Convergence of different fidelity algorithms in a 100- 
dimensional system of perturbed (e = 3 x 10~ 4 ) quasi- 
integrable (k — 0.2) kicked rotors with m = 8192. Error bars 
plotted every 20 time steps, (a) Simple algorithms echo-1 and echo- 
1' are far from converged even with 7 x 10 7 trajectories, (b) While 
both DR and echo-2 algorithms converge fully with only 2048 trajec- 
tories, only the DR can capture the QM fidelity "freeze" (the plateau). 



FIG. 2. Statistical error grows exponentially with D for the echo-1, 
echo-1', and echo-N-1 algorithms, while it is independent of D for 
the echo-2 and DR algorithms, (a) Pure displacement dynamics ob- 
tained with two displaced D-dimensional SHOs. TV ~ 10 7 . Time 
chosen separately for each D so that F ~ 0.3. (b) General dynam- 
ics obtained with a .D-dimensional system of perturbed (e = 10~ 4 ) 
quasi-integrable (k — 0.2) kicked rotors with n\ = 131072. N ~ 
5 x 10 5 . Time chosen separately for each D so that F « 0.9. 



for <7^ dr (exact) and u Fbr (to leading order in N) gives 



N fDR =a- 2 (l-F DR )and 



Nf = — 



Re 



(<■ 



i24>\ 



I Pw ^ 



-i<j>\ 



2F" 



(14) 



(15) 



Result (O for N 



fm 



is completely general. As for Np DR , using 
< 1 and Eq. (O, we can find a 
completely general upper bound, 



the inequality | (e 42<#, ) pw(x o ) 



Nf dr < 4<J- 2 F DR (l - Fi 



DR) 



(16) 



Estimate ( TBI and upper bound (fl~6] i show, remarkably, that 
without any assumptions, the numbers of trajectories needed 
for convergence of both Jur and Fur depend only on a and 
F DR , and are independent of D, t, initial state, or dynamics. 
Estimate ( fl~5b of Np DR can be evaluated analytically for nor- 
mally distributed phase tj>. This is satisfied very accurately in 
the chaotic FGR and integrable Gaussian regimes , and 
exactly for pure displacement dynamics of GWPs. Noting that 
for normal distributions (e 4 *) = e l '^'exp[— Var((^>)/2] and 
Fm = | /dr = cxp[— Var(0)], Eq. ( fl3] > reduces to 



N, 



2a- 2 F DR (1 



Fdr) ' 



(17) 



which is again independent of D, t, initial state, or dynamics. 



Using a similar analysis, in Ref. [4] it was found 
that for CF algorithms ([T0li-(fT3]> and D 1, one needs 
N = a~ 2 a(F)f3 D trajectories where a and j3 depend on the 
method, initial state, and dynamics. For all methods with 
M 7^ 2, there are simple examples 10] with j3 > 1, imply- 
ing an exponential growth of N with D. Remarkably, for any 
dynamics and any initial state, (3 = 1 for the echo-2 algorithm, 
implying, as for the DR, that N is independent of D [4]. 

Numerical results and conclusion. To illustrate the analyt- 
ical results, numerical tests were performed in multidimen- 
sional systems of uncoupled displaced simple harmonic os- 
cillators (SHOs, for pure displacement dynamics) and per- 
turbed kicked rotors (for nonlinear integrable and chaotic dy- 
namics). The last model is defined, mod(27r), by the map 
Qj+i = <h ■ l' r /', • i = Pi -VW(<fr+i)-eVV(g i+1 ) where 
W(q) = —kcosq is the potential and V(q) = — cos(2<?) 
the perturbation of the system; k and e determine the type of 
dynamics and perturbation strength, respectively. Uncoupled 
systems were used in order to make QF calculations feasi- 
ble (as a product of D 1-dimensional calculations); however, 
both CF and DR calculations were performed as for a truly 
D-dimensional system. The initial state was always a multi- 
dimensional GWR Expected statistical errors were estimated 
by averaging actual statistical errors over 100 different sets 
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FIG. 3. Regardless of dynamics, statistical error of the DR is in- 
dependent of dimensionality (D), time (t), and is proportional to 
iV -1 / 2 . Errors are compared for 10 kicked rotors in the chaotic FGR 
regime (k = 18, e = 6.4 x 1CT 6 , m = 131072), 100 kicked rotors 
in the integrable Gaussian regime (k — 0.2, e = 6.4 x 10 _6 ,ni = 
131072), and a single kicked rotor in the quasi-integrable algebraic 
regime (k = 0.2, e = 6.4 x 10" 4 ,ni = 131072). Timet was 
chosen separately for each system so that F ~ 0.94. 



of N trajectories. No fitting was used in any of the figures, 
yet all numerical results agree with the analytical estimates. 
Note that the figures show also results for algorithm echo-1', 
■Fecho-r (t) = 1 + {p{x~ l ) - p(x°)h^ x oy which is a variant of 
echo- 1 accurate for high fidelity Jj] . 

Figure Q] displays fidelity in a 100-dimensional system of 
kicked rotors. It shows that both echo-2 and DR algorithms 
converge with several orders of magnitude fewer trajectories 
than the echo-1, echo-1', and echo-N-1 algorithms but while 
the DR agrees with the QM result, even the fully converged 
CF (computed as a product of 100 one-dimensional fidelities) 
cannot reproduce QM effects. Figurel2lco nfirms that whereas 
the statistical errors of the echo-1, echo-1 , and echo-N-1 al- 
gorithms grow exponentially with D, stat stical errors of the 
DR and echo-2 algorithms are independent of D. Figure [3] 
shows that for several very different dynaj nical regimes, ctdr 
is independent of t, D and n\, in agreement with the general 
upper bound ( [ToT i and-in the FGR and Gau ssian regimes-also 
in agreement with the analytical estimate 



ureg]exhibits the superior computational efficiency of the DR [3] R Jac <3 uod and c - Petitjean, |A.dv. Phys., 58, 67 (2009 

!V|lollica, T. Zimmermann, 
[108.0173 [nlin.CD] 



compared to all CF algorithms: thanks tc 



with t and independence of D, the DR is o rders of magnitude 
faster already for quite a small system and short time. 

To conclude, in the case of QF, a SC nethod can be not 
only more accurate, but also more efficient than a CL simu- 
lation of QD. This counterintuitive result s lould be useful for 
future development of approximate metho ds for QD of large 
systems. This research was supported by S viss NSF grant No. 
200021.124936 and NCCR MUST, and by EPFL. We thank 
T. Seligman and T. Zimmermann for discussions. 
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